Stockname = char('AAA');
Indexname = char('HNINDEX');
[Stock Index Date PriceStock PriceIndex] = getStock_Index('AAA','HNXINDEX');
%Get the return of stock and index
rStock = PriceStock(2:end)./PriceStock(1:end-1) - 1;
rIndex = PriceIndex(2:end)./PriceIndex(1:end-1) - 1;
%-----------------------------------------------------------
%---Estimator: OLS 1Year-BETA BAYESIAN
%-----------------------------------------------------------
%Rolling window
rolling = 260;
[n,m] = size(rStock);
k = n-rolling;
datan = x2mdate(Date,0);
datam=datan(262:end);
    InrStockOls=rStock(1:260,:);
    InrIndexOls=rIndex(1:260,:); 
    [n,m] = size(InrStockOls);
    for i=1:m
      bo(i) = regress(InrStockOls(:,i),InrIndexOls);
    end 
    betapr = mean(bo);
    varbetpr = std(bo)



for i=1:k
    j=i+260;
    rStockBayes = rStock(i:j,1);
    rIndexBayes = rIndex(i:j,:);
    [n,m]=size(rStockBayes);
    X= [ones(n,1),rIndexBayes];
    y=rStockBayes;
    %QR Decomposition X
    [Q,R]=qr(X,0);
    %Regression Coefficients 
    beta = R\(Q'*y);
    %Fitted Values of the Response 
    yhat = X*beta;
    %Residuals
    residuals = y - yhat;
    residualsm(i)=mean(residuals);
    %Mean Squared Error 
    nobs = length(y);
    p = min(size(R));
    dfe = nobs-p;
    mse = sum(residuals.*residuals)./dfe;
    sse(i)=sum(residuals.*residuals);
    mse(i)=mse;
    %Hat (Projection) Matrix 
    hatmat = Q*Q';
    yhat = hatmat*y;
    %Leverage 
    leverage = diag(hatmat);
    %Covariance Matrix of Estimated Coefficients 
    ri = R\eye(p); % inverse of R
    xtxi = ri*ri'; % equivalent to inv(X'*X)
    covb = xtxi*mse(i);
    %Student's t statistics 
    coeff(i,:) = (R\(Q'*y))';
    se(i,:)= sqrt(diag(covb));
    t(i,:)= coeff(i,:)./se(i,:);
    dfe=nobs-p;
    pval(i,:)= (tcdf(-abs(t(i,:)),dfe)).*2;
    %F statistic 
    sse(i)  = norm(residuals).^2  ;
    ssr(i)  = norm( yhat - mean( yhat)).^2;
    fdfe  =  nobs-p;
    dfr  = p-1;
    f(i) = (ssr(i)/dfr)./(sse(i)/dfe);
    pfval(i) = 1 - fcdf(f(i),dfr,dfe);
    %Stima Bayes
    meanind=mean(rIndexBayes);
    varres(i)=var(residuals);
    v(i)=sum((rIndexBayes(:,1)-meanind).^2);
    varols(i)=varres(i)/v(i);
    varbayes(i)=1/(1/varols(i)+1/varbetpr);
    w(i)=varbayes(i)/varbetpr
    betabayes(i)=w(i)*betapr+(1-w(i))*coeff(i,2);
end

%Output to Excel
ExcelBayes=[datam,betabayes'];